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Abstract. We study an analytically tractable model with long-range interactions for which an out-of- 
equilibrium very long-lived coherent structure spontaneously appears. The dynamics of this model is indeed 
very peculiar: a bicluster forms at low energy and is stable for very long time, contrary to statistical me- 
chanics predictions. We first explain the onset of the structure, by approximating the short time dynamics 
with a forced Burgers equation. The emergence of the bicluster is the signature of the shock waves present 
in the associated hydrodynamical equations. The striking quantitative agreement with the dynamics of the 
particles fully confirms this procedure. We then show that a very fast timescale can be singled out from a 
slower motion. This enables us to use an adiabatic approximation to derive an efi'ective Hamiltonian that 
describes very well the long time dynamics. We then get an explanation of the very long time stability of 
the bicluster: this out-of-equilibrium state corresponds to a statistical equilibrium of an effective mean-field 
dynamics. 
Keywords: 
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1 Introduction 

The Hamiltonian Mean Field model (HMF) has attracted 
much attention in the recent years as a toy model to study 
the dynamics of systems with long-ranqe interactions, and 
its relation to thermodynamics The HMF model 

describes an assembly of N fully coupled rotators, whose 
Hamiltonian is: 



N 



2N 



N 

E 

2J = 1 



cos 



(1) 



where 9i is the angle of the i-th planar rotator with a 
fixed axis. As the interaction only depends on the angles 
of the rotators, this model can alternatively be viewed as 
representing particles that move on a circle, whose posi- 
tions are given by the 9i and interact via an infinite-range 
force. If c is negative, the interaction among rotators is 
ferromagnetic, corresponding to an attractive interaction 
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among particles. When c is positive, the interaction among 
rotators is antiferromagnetic and repulsive in the particle 
interpretation. In this article, we will focus our study on 
this latter case, and put c = 1. 

As first noticed in Ref. this model has a very in- 
teresting dynamical behavior. In contrast to the statistical 
mechanics predictions, a bicluster forms at low energy (see 
Fig. ^, for a special, but wide class of initial conditions: 
it appears as soon as the initial velocity dispersion of the 
particles is small, for any initial spatial distribution. This 
clustering phenomenon and its unexpected dependence on 
initial conditions were precisely studied in Ref. @|, but 
remained essentially unexplained. Several details of the 
dynamics were numerically studied and it has been in 
particular emphasized that the energy temperature re- 
lation is modified: although still linear, the slope in the 
molecular dynamics simulations is different from the the- 
oretical prediction. No sign in the numerics was found of 
the decay to the constant density profile predicted by the 
canonical ensemble. On the contrary, recently performed 
simulations |^ with a smaller number of particles have 
shown a long-time degradation of the bicluster, suggest- 
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ing its transient non-equilibrium nature. The question of 
the time asymptotic stability of the bicluster, in the limit 
of an infinite number of particles, remains however open 
and will be discussed in the conclusions. 

We will consider in this work two different theoreti- 
cal approaches that lead to an analytical explanation of 
the clustering process and of the long time stability of the 
structure. In section in order to explain the short time 
bicluster formation, we first propose a hydrodynamical de- 
scription that links Hamiltonian ([l]) to a forced Burgers 
equation (these results were shortly presented in a pre- 
vious brief note This striking phenomenon is then 
described in Lagrangian coordinates, giving rise to inter- 
sections of trajectories and formation of high (or even in- 
finite) particle densities. We then solve in section |^ this 
equation, using extensively the method of characteristics. 
Beyond the theoretical interest in such structures for Hamil- 
tonian systems, we will show some connections of this 
problem with other subjects; namely, active transport in 
hydrodynamics, and the formation of caustics in the flow 
of the associated Burgers equation. 

Section ^ presents the second approach, which is par- 
ticularly powerful for explaining the long-time evolution 
of the bicluster. The dynamics involves two well sepa- 
rated time-scales. Using a variational approach, we will 
show how to construct an effective Hamiltonian where all 
the fast oscillations are averaged out. Not only this effec- 
tive dynamics accurately represents the one of the origi- 
nal Hamiltonian but, in addition, the complete statistical 
thermodynamics can be easily derived in the microcanon- 
ical and canonical ensembles. It predicts the existence of 
the bicluster as an equilibrium state, with the correct den- 
sity profile and the microcanonical temperature energy 
relation found in the numerical experiments of Ref. B. 
The long lifetime of these out-of-equilibrium states can 
therefore be interpreted by the fact that they appear as 
equilibrium states of an effective Hamiltonian represent- 
ing the long-time motion. Section || is devoted to some 
conclusions and perspectives of future developments. 
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Fig. 1. Bicluster formation : short-time evolution of the par- 
ticle density in grey scale: the darker the grey, the higher the 
density. Starting from an initial condition with all the parti- 
cles evenly distributed on the circle, one observes a very rapid 
concentration of particles, followed by the quasi periodic ap- 
pearance of "chevrons", that shrink as time increases. 



and a velocity field 

Pi9,t)v{e,t) 



+ 00 



pfi9,p,t)dp 



(4) 



2 Hydrodynamical description at zero 
temperature 

2.1 Introduction 

To describe the initial clustering process, we consider the 
associated Vlasov equation, which can be rigorously de- 
rived in the thermodynamic limit N — > oo Denot- 
ing by f{9,p,t) the one particle distribution function, we 
have: 



a/ 5/ 
at ^86 
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2^ 
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du / da /(a, u, i) sin(6' — a) 



Let us define as follows, a density 
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fie,P,t)dp 



di 

dp 
(2) 



(3) 



As the numerical simulations reported in Ref. have 
shown that the bicluster appears when the velocity disper- 
sion is small, we will consider here the zero temperature 
approximation, i.e. we neglect the velocity dispersion and 
consider the ansatz f{6,p,t) = p{9,t)S{p — v{9,t)), where 
S refers to the Dirac function. The Vlasov equation (||) 
can consequently be simplified to 



5/ 
dt 
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where one notes 
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The sum of the time derivative of Eq. 
derivative of Eq. (m , leads to 



(5) 



(6) 



and the 9- 



dp djpv) 
dt 39 
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(7) 



We obtain therefore the equation which accounts for mass 
conservation 

dp djpv) _ 

Using the time derivative of Eq. (||) and Eq. (H), we 
obtain 



d{pv) 



dv dv 



d{pv'^) 

de 



(9) 



However, the time derivative of the left-hand-side gives 
also 



dt 



(10) 



Equations and Eq. ( [lo| ) lead finally to the following 
Euler equation without pressure term: 



^^""^6 = 2^1 d-PMMO-a) (11) 

Equations and ( pi] ) show that the dynamics of the 
model at low temperature can be therefore mapped onto 
an active scalar advection problem. 



2.2 Linear analysis 

Linearizing equations (||) and ( pl) , by assuming small ve- 
locities (f <C 1) and almost uniform density (p = 1 + 
pi{0,t) with pi <^ 1), we are left with: 



dv 



de 



dpi 
dt 

dv_J_ 
dt ~ 27r 



= 



2ir 



Pi (a, t) sm{0 — a) 



(12) 
(13) 



a system which is easily solved by a spatial Fourier series 
development. Assuming a zero initial velocity field, we get 

Pl{9,t)=^/2 VmCOsd COBUOt (14) 

t;(6', t) = u,„ sin6' sinwt , (15) 

where the time-scale of the oscillation is given by the in- 
verse of a sort of "plasma frequency" lo = •\/2/2, which 
describes, as in the Poisson-Boltzmann case, the instabil- 
ity of the uniform density state . In the attractive case 
(c = —1 in Hamiltonian (|l|)), this time scale would con- 
trol the depart from the initial uniform density towards 
the formation of the density profile of the single cluster 
that appears at low temperature [ pO| , 1 1 1 , in analogy with 
Jeans instability in gravitational systems [ p^ . 

This linear analysis is, however, not sufficient to ex- 
plain the formation of the bicluster and we have to carry 
out a non linear analysis. 



2.3 Non linear analysis 

This analysis relies on the existence of two time-scales in 
the system: the first one is intrinsic and corresponds to 
the inverse of the "plasma frequency" the second one is 
connected to the energy per particle e. When the energy 
is sufficiently small, these time-scales are very different, 
and it becomes possible to use averaging methods. The 
presence of two well separated time-scales also explains, 
in some sense, the appearance of the bicluster in the low 
energy limit. 

We introduce therefore a long time-scale r = et, with 
£ — Vm/V^ — \/2e, and we look for solutions of the fol- 
lowing form: 



v{e,t) 



, sin 6 smujt - 



uie,T) 



(16) 



In order to evaluate the force on the r.h.s. of Eq. (|11|) in the 
nonlinear regime, we will use expressions (y) and (pj), 
given by the linear analysis. Indeed, due to the special 
form of the Hamiltonian, the force depends only on the 
first Fourier component of the density. Hence, our hypoth- 
esis amounts to assume that the sinusoidal behavior of the 
density, found in the linear regime, holds the same in the 
non-linear regime: the results presented below confirm the 
validity of this assumption. 

In analogy with studies of the wave-particle interac- 
tion in plasma physics |l^,|lj|, our system may be seen 
as a bulk of particles interacting with waves that are sus- 
tained by the bulk itself. Here, the wave is created by the 
small density and velocity oscillations already present in 
the linear analysis. In the low energy regime, the phase 
velocity of the wave with the plasma frequency uj is much 
higher than the velocities of the bulk particles, causing 
a very small wave-particle interaction strength: the wave 
has therefore a very long life-time. This also explains the 
quality of the linear approximation for the force that we 
have used in the Euler equation (pT|). 

Introducing expression ( |l6| ) in Eq. (|ll]), terms to first 
order in e disappear by construction. Order terms give: 

du du „ . „ /, . 2 

t: h u^rr + 2smt^coswsm uit 

dr dO 

sin 6*^ +MCOS0 I \/2sincjt = (17) 
do J 

Averaging over the short time scale t, we obtain: 



du du 1 . „^ 

+ = 



(18) 



which is a spatially forced Burgers equation without vis- 
cosity, describing the motion of fluid particles in the po- 
tential V{e) = 1/4 cos 26*. 

The Burgers equation has been studied in many dif- 
ferent contexts by mathematicians interested in a pres- 
sureless description of fluid motion, and by physicists in 
the context of ballistic aggregation models . In partic- 
ular, the forced Burgers equation has been very carefully 
studied from a mathematical viewpoint |l^ . Cosmologists 
have proposed a very interesting adhesion model that uses 
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Burgers-like equations to explain structure formation in 
the universe Moreover, this equation may be also 

viewed as a nice model of active transport of the vorticity 
field in fluid mechanics fl^ . 

A well known property of the Burgers equation with- 
out viscosity, is that the solution becomes multi-stream 
after a finite time: the appearance of shocks (see Fig. ||) 
in the velocity profile u(9) corresponds to the creation of 
singularities in the density profile (see Fig. ||). This is a 
consequence of the two main assumptions: the medium is 
supposed to be continuous and the temperature is set at 
zero. Weakening either assumption would have smoothed 
out the singularities. In the original discrete Hamiltonian 
model, particles can cross and, after some time, those that 
travel faster can be catched by the slower ones that lie 
downstream, creating the spiral dynamics exemplified in 
Fig. (^. The particle density has a peak at the center 
of the spiral (Fig. ||), and diverges at the peak position 
in the N —^ oo limit at the shock time r^. The presence 
of a shock at finite time in the forced Burgers equation 
does not prevent the hydrodynamical description from re- 
maining valid at longer times, when, moreover, several re- 
peated shocks are observed at regular time intervals. We 
will indeed show that it can adequately describe even the 
long-time evolution. In addition, let us remark that it's the 
double well shape of the potential that forces the Burg- 
ers equation, which is responsible for the formation of two 
clusters. Starting from the initially uniform state, the par- 
ticles will move around the bottom of the both wells and 
the majority of them happen to be there at the same time 
at the shock time, creating therefore two coherent struc- 
tures that after several oscillations form the bicluster at 
long-time. 
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Fig. 2. Shock Dynamics. Phase space portrait of = 10* 
particles, which at r = \/2ei ~ are uniformly distributed 
in space (dotted line) with a small sinusoidal velocity profile, 
not visible in the figure. The resulting energy is e = 6.7 10"'*. 
Velocity profiles u{9) are then shown at r = 7r/4 (dashed line), 
r = 3tt/8 (dashed-dotted line) and t — Ts — 7r/2, the first 
shock time (solid line) . Only half the space is shown since the 
curves are 7r-periodic. 



3 Solution of the spatially forced Burgers 
equation by the method of characteristics 

The method of characteristics is the appropriate mathe- 
matical tool to study more quantitatively the forced Burg- 
ers equation ( p^ ) , but it has also a direct physical meaning: 
the characteristics correspond to the Lagrangian trajec- 
tories of the Euler equation, and are therefore good ap- 
proximations for the particle trajectories of the finite N 
Hamiltonian system (|^) when TV ^ 1. 

In the case of unidirectional nonlinear wave motions, 
the method is standard and proceeds as indicated, for ex- 
ample, in Ref. We obtain in this case the following 
pendulum equation for the characteristics 0{t): 

+ isin2^ = O, (19) 

with the initial conditions 6(0) = 9o and ^(0) = 0. The 
solution of the above equation can be easily derived |^ 
and reads: 

= arcsin [sin 00 sn (r + Xq, sin^o)] (20) 

where — K{sm9o), K is the complete elliptic integral 
of the 1st kind and sn, the elliptic sine function. Fig. (^) 
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Fig. 3. Particle density p{9) at r = (dotted line), r = 7r/4 
(dashed line), r = Stt/S (dashed-dotted line) and r = Ts = 7r/2 
(solid line). The initial condition is the same as in Fig. ^ We 
clearly see the appearance of one of the two density peaks of 
the bicluster. 
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Fig. 4. Spiral Dynamics. Phase space portrait of A'' = 10* 
particles at r = 500. The initial condition is the same as in 
Fig-i 

shows the characteristics, i.e. the solutions of Eq. ( p^ 
corresponding to different initial conditions evenly dis- 
tributed in [0, 2tt\ at t = and with zero initial velocity. 
One clearly sees their oscillatory behavior inside one of the 
two wells of the effective potential cos 20, located at = 0, 
with a period 'iK^ which strongly depends on the initial 
position. Fig. (|^) emphasizes also that the characteristics 
cross themselves, and the associated caustics are respon- 
sible for the enhancement of the particle density field that 
forms the "chevrons" . The time of the first divergence in 
the density Tg is the time of the first characteristic crossing 
and also the time of the first shock in the velocity profile. 
The particles performing a quasi-harmonic oscillatory mo- 
tion in the bottom of the well are the ones with shorter 
periods. The shortest period is obtained in the harmonic 
limit 4i^(0) = 27r, leading to Ts = 7r/2 for the time of ap- 
pearance of the first shock. In the original units, it reads 

ts-"-Z-^ ■ (21) 
£ V8e 

This shock results in a singularity in Eulerian space 
at the corresponding time. However this singularity dis- 
appears immediately after forming and two singularities 
of another kind arise in its place, which are the bound- 
aries of the three streams region. The chevron is their 
manifestation in the density profile. The shock time sin- 
gularity is called in Arnold's classification |^,^ while 
the "chevrons" singularity is of A2 type. The latter is the 
boundary of structures and can exist at any moment of 
time, whereas the former exist only when a chevron origi- 
nates. The recurrence time for the appearance of chevrons 




2 4 6 8 10 



T 




Fig. 5. Panel (a) presents the characteristics of equation ( |l8| ) 
whereas panel (b) presents the trajectories of the particles of 
the Hamiltonian (Q) with the same initial condition as in Fig. ^ 
One can see the three first appearances of the "chevrons" . Two 
phenomena are not captured by the characteristics: the small 
oscillations of the real trajectories, which are averaged out, and 
the presence of untrapped particles, close to the saddle-points 
of the effective cos 26 potential. 
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is t„ = ts * (2n— 1) and the comparison with the numerics 
shows a very good agreement. 

Figs. (I^a) and (||b) show that the agreement is reaUy 
good not only at the qualitative but also at the quantita- 
tive level. Only two features are missed by the Lagrangian 
description. The fast oscillations of very small amplitude, 
already hardly visible in Fig. (||b) is of course totally ab- 
sent in the characteristics because of the averaging proce- 
dure used to get Eq. ([l8|). One can also easily show that 
the amplitude of this fast oscillation is rapidly decreas- 
ing with the increase of the number of particles and is 
vanishingly small for N oo. The second point is the 
presence of untrapped particles in Fig. ^d), which is a 
direct consequence of the oscillations of the height of the 
potential barriers. This effect is important only for highly 
energetic particles whose trajectories in phase space are 
closed to the separatrix of the Kelvin's cat eye of Fig. . 
However, as we will show in the next section, these par- 
ticles do not take part to the creation of caustics (i.e. of 
the chevrons of Fig. |l|) which are generated only by the 
particles close to the bottom of the wells. 



3.1 The chevrons as caustics of the characteristics 

We will now explain the shape of the " chevrons" (Fig. |l|) , 
which correspond to zones of infinite density i.e. to the 
envelops of the characteristics, the so-called caustics. The 
family of the characteristics is defined by 



— sin 9o and we obtain thus 



F{T,d, 



%) - sin(6l) = 



(22) 



in the plane {O^t) with the parameter 9q. The envelope 
of this family would then be defined by the two following 
equations 



F(t,0,0o) = 
o) = 



dF , 



(23) 
(24) 



However, it is not possible to extract 9q as a function of 9 
and T from Eq. ( p^ ) , in order to obtain a closed expression 
F{t,9,9o{t,9)) for the caustics. We will use approximate 
expressions of F close to the shocks. 

In the neighborhood of the shock, the trajectories can 
be approximated by straight lines, with the following ex- 
pression: 



(r, 6*0) = Vq (t - Kq 



(25) 



where vq = v{Kq) is the speed of the particle at the bot- 
tom of the potential well. 

The envelope of this family of characteristics can then 
be obtained Q by solving the following system of equa- 
tions 



(r,0o) (t-Ko) 

= ^'o (t - Kq) - vqKq 



(26) 
(27) 



where the primes denotes the derivative with respect to 
Oq. Using the conservation of energy, we easily get vq — 



T — tan 9o Kq + Kq 
_ sin^ 9q , 



(28) 
(29) 



It is however possible to go further by using higher 
order terms in the expression of 9 given in Eq. (p5|), i.e. 
by considering curves rather than simply straight lines. 
We have thus 



(r,9Q) = 9{KQ) + ^ (t-Ko) 

dT\Ko 

1 S9 



6 dr^ \Ka 

VQ (r-KQ)-^iT-KQf 
6 



2 dr"^ \Ko 

- Kq)' -f 



(r - Ko)' 



(30) 
(31) 



where we have replaced in Eq. ( pi] ) the derivatives by their 
expressions. Let us note in particular, that all derivatives 
of even order will vanish because of the parity property 
of 9. Using the development of Kq at the same order, 
we finally end with the next order approximation for the 
caustics 



3a9l 



5b 



7 c 



9 = 2a9, 



Ab 

y 

2b 



2a 

y, 

4a 8a^ 
15 ^ ~3~ 
a 

3 

31a , , 
180 



0{9l) (32) 



+ 6c+^ + ^-4a3 0^ + 0(^9) (33) 



where (a, 6, c) 



11 173 
48' 2880 



It is possible to see in Fig. (g) that the agreement is 
excellent and that this procedure is really accurate to de- 
scribe the first chevron. Moreover, one can check that be- 
cause of the non isochronism of the oscillations, the more 
energetic particles arrive too late in the bottom of the well 
to take part to the creation of the caustics, i.e. the singu- 
larity in the density space. Consequently, the untrapped 
particles of the real dynamics, visible in Fig. (||b), but ab- 
sent of the averaged dynamics of the associated Burgers 
equation (see Fig. |^) , are irrelevant and do not affect the 
dynamics of the caustic formation. 

This is confirmed by Fig. (^ where the caustics de- 
termined by the Burgers' approach are superimposed on 
the real trajectories: the shape of the chevrons is very well 
described. However, the time occurrence of the shock was 
reduced by a factor 0.84. The reason, clarified in the next 
section, is due to the renormalization of the oscillation fre- 
quency uj which modifies the oscillation period of the par- 
ticles. As shown below using a self-consistent procedure, u) 
is reduced with respect to the plasma frequency ujp. As the 
potential height is inversely proportional to uj this leads 
therefore to a faster appearance of the shocks. 

As already mentioned, the time occurrence of the nth 
shock is T„(0o) = (2?T- — 1)ts which means that its shape 
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Fig. 6. Superposition of the caustics over the characteristics of 
50 particles evenly distributed between —n/2 and n/2 att = 0. 
The analytical formula for the "chevrons" (Eq. ^ and Eq. ^3|) 
is superimposed (bold curves) . The vertical dashed lines shows 
the appearance time of the chevrons tn = ts * (2n — 1). 
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Fig. 7. Superposition of the caustics over the trajectories of 
100 particles evenly distributed between and 27r at t = for 
an energy e = 3 10~®. The time occurrence of the shocks was 
reduced by a factor 0.84 (see text). 



will be obtained simply by replacing in Eq. ( |2^ ) Kq by 
(2n — 1)Kq. Therefore, the lowest order approximation of 
the nth shocks is 



9 = 2a(2n- 1) 



,3a(2n - 1) 

(^-(2n^l)t,)^/' 
V2n- 1 



3/2 



(34) 
(35) 



where n is the number of appearance of the "chevron". 
The 1 / y/n factor accounts for the shrinking of the "che- 
vrons" and the bolded line in Fig. (||) attests that this 
expression is particularly accurate. 

Similar caustics are encountered in astrophysics, to ex- 
plain the large scale structure of the universe: clusters and 
super clusters of galaxies are believed to be reminiscent of 
three dimensional caustics arising from the evolution of 
an initially slightly inhomogeneous plasma p3,E3]. 



3.2 Probability distribution 

At this point, it would be important to calculate the den- 
sity distribution in the vicinity of the singularity. For this 
hydrodynamical approximation, it would be fully deter- 
mined by the initial conditions, i.e. the initial velocity dis- 
tribution. Let us recall that we earlier found [Q, that the 
following analytical formula 



1 

2^ 



(l-log(2|sin0|)) 



(36) 



was an excellent approximation (see Fig. even if we 
didn't have any theory for its derivation. Here using the 
Burgers' approach, let us obtain a similar result. 
Along the particle trajectory, one has 



cos 29 — cos 200 



(37) 



The time dt{9) spent by a particle close to a position 9 is 
inversely proportional to its velocity 9{9q). Therefore the 
trajectory parametrized by 9q (and initiated in ^o) gives 
a contribution to the density 



PeMde 



dt{9)d9 
dt{9) 



2^/2 d9 



T{9o)y/cos29 -cos29o 



(38) 



where T{9o) is the period of the trajectory. By recaUing 
that the initial distribution is homogeneous in these nu- 
merical computations, one obtains the probability density 
by averaging over the time, the density of characteristics 
at a fixed position. We obtain: 



p(0)(x 



2^/2 du 



28 



T{u)y/cos29 - cos2u 



(39) 



for 9 = [0, 7r/2] whereas the whole distribution is obtained 
by symmetry and 7r-periodicity. 
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Fig. 8. Equilibrium distribution. Comparison between the 
numerical results (diamonds) and the analytical formula ( p9[ ) 
(solid line) in the case e ~ 10~*. The dotted line corresponds 
to the formula (pd). 



Due to the numerous approximations, the agreement 
for long time is not perfect, but equation (|9|) gives a good 
result as shown by Fig. (H) . The disagreement is visible in 
and TT, i.e. close to the separatrix. Actually, as we will see 
in the following, the amplitude of the effective potential 
slightly oscillates; this allows the existence of untrapped 
particles which smooth the density. 

As a conclusion, this hydrodynamical approach has 
proved to be very successful, since it explains qualitatively 
and quantitatively, for short times, the main features of 
the bicluster formation. Nevertheless, the stability and the 
long time behavior of the structure is not yet understood: 
this is the object of the next section. 



4.1 The fast variables 

The complete Lagrangian of the antiferromagnetic HMF 
being 



^('.■^O-St-^^i:' 



(40) 



the equations of motion are obtained by minimizing the 
action S = J Ldt with respect to the functions 9i. Taking 
advantage of the two well defined time scales, we introduce 
the following ansatz: 



9Ur)+sMt, 



(41) 



wher e r — et and e = v2e as already defined in the sec- 
tion (pT^). The fi's represent the small rapid oscillations 
and the 0^ the slow variables which we are finally inter- 
ested in. The idea is then to use a variational approach 
(least action principle), first to obtain the equations of 
motion of the fast variables fi, and then to reintroduce 
the solutions into the action. In a second stage, averaging 
over the fast time t, we will obtain a variational princi- 
ple for the slow variables. This approach, inspired by [p4| , 
has the double advantage of allowing to obtain an Hamil- 
tonian system for the slow variables, and the one of ex- 
hibiting a natural conservation law for this dynamics (the 
adiabatic invariant). Appendix A briefly presents the use 
of method for a slowly modulated harmonic oscillator, 
in order to emphasize its power for a very simple model, 
away from the rather technical context presented here. 

We first notice that since the energy of the system 
is by definition of order the sums l/iV^cos^^ and 
1/A^^sin0,°, which represent the two components of the 
magnetization vector in terms of the slow variables 0^, are 
of order e. We thus define 



i 



(42) 
(43) 



4 Long time evolution of the bicluster 



where the phase ip G [0, tt] is chosen such that the scalar 
dynamical indicator of the clustering is 



Since the thermodynamics of the antiferromagnetic HMF 
model predicts a perfectly homogeneous equilibrium state, 
the bicluster is a non equilibrium structure. However, its 
long time stability, as well as the fact that it is reached 
starting from a variety of initial conditions suggest some 
statistical explanation. 

In addition, we have learned from the previous sec- 
tion that particle motion may be split into two parts: a 
small amplitude rapid oscillation superimposed to a large 
amplitude slow motion. The idea is thus to look for an 
effective Hamiltonian dynamics which would only retain 
the slow motion, and would be appropriate for a statistical 
description. 



|M2 



l^exp(2^0°) .l^cos2(0O+^). (44) 



We now introduce the ansatz (^Tj) into the Lagrangian (pO|), 
and develop the cosine up to order e^, obtaining the new 
Lagrangian L2 that depends on the the fi's and their 
time derivatives 



L-2 — L2 



del 

dr 

del 

dr 



fi-! 



dt 



dr dt 
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Due to the averaging, the variables (j)± are cyclic, P±, the 
conjugate momenta of (f>±, are conserved quantities. Their 
expression is 



(45) 



(54) 



From this expression, considering r as a constant, we write 
the equations of motion for the fast variables /^'s. We get 



/, = Ml sin (00 + V-) -Ml cos (0° + ^p) 



As the Lagrangian does not depend on the time derivatives 
of , A- , a+ , and a_ , there is no Legendre transform 
on these variables and the Hamiltonian reads 



; cos I 



(46) 



N 



Heff = P+ 



This is a linear second order equation, with constant coef- 
ficients with respect to the fast time, whose solution thus 
requires the diagonalization of the NxN matrix 



Pi 



p2 



30 



N q2 
Pi 



(55) 



[rd = ^[cos(0°-^°)] , (47) 
which has only two non zero eigenvalues (see Appendix B) 

2 1±|M2| 



2NAl'''^ 
+N 



2 



E 



Ml 



+ Mla+uj^ - Mla^ujt 



(48) 



2 



A' 



(56) 



The eigenvectors corresponding to are respectively 



X+ = [cos(0O + V^)]^^^^ ^ 
X_ = [sin(0,? + V')],^i...^ . 
The general solution for the fiS is therefore 

il 



(49) 
(50) 



In the absence of conjugate variables of the amplitudes 
A±, the corresponding Hamilton's equation are simply 
given, from the least action principle, by dA±Heff — 0. 
Together with the equations (^4|) , this leads to the follow- 
ing expressions for the frequencies 



V2A+ sin(a;+t + (/?+) — 
V2A^ sin(w_i + (^_) + 



Ml 



sin(6i° 



V;) 



dt 



(57) 



(51) 



expected of course from the previous study of the ma- 
trix T, and for the amplitudes 



4.2 The slow variables 

The previous solution for the fi suggests the following new 
ansatz 



/, = [V2A+(T)sin((^+(i)) + a+(r) 
+ [\/2A_ (t) sin(0_ (t)) +a_(r) 



cos 



9." + V') 
?^+^) (52) 



Finally, from the equations da± H^f f 
M^ 

a+ = if- and a_ = 



= 0, we find 

Ml 



(58) 



(59) 



where 4>±{t) are fast variables, and d4)±/dt, A±{t), and 
a±(T) are slow variables. We introduce ( p^ ) in the La- 
grangian (^), and we average over the fast variables 0±. 
Dropping the overall factor, the averaged Lagrangian 
reads 



Reintroducing expressions (^8|) and (^) in the Hamilto- 
nian (p6[), we end up with the effective Hamiltonian de- 
scribing the slow motion of the particles 



ff 



P^UJ^ 



(60) 



\ rl^ ^9 



N 



i=l 



-N 



N 



At 



MV + Ml 



u\ - Mla^ 

A 



+ {Al + al)'^ + {Al+al) 



,(53) 



The evolution of the full system is then approximated 
by the dynamics of this effective Hamiltonian, the con- 
stant P+ and P_ being determined by the initial condi- 
tions. However, a difficulty arises when the two eigenvalues 
crosses, ie when IM2I becomes 0. Since -0 is defined as half 
the phase of the complex number M2X + iM2y, it expe- 
riences at this point a 7r/2 jump, and consequently the 
eigenvectors X± are inverted: this phenomenon of eigen- 
value crossing is illustrated by Fig. ^. 
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Fig. 9. The diamonds presents the evolution versus time of the 
main frequency for — 200 particles of the original Hamilto- 
nian (0). We have superimposed not only uj+ (dashed line) 
and Lj^ (solid line) using Eqs. (^^, but also the numerically 
computed IM2I (dotted line). We clearly see the exchange of 
frequency when \M2\ touches zero. 



However, in most of the numerical experiments con- 
ducted in B, only one fast frequency, uj- was excited. 
From now on, we will therefore concentrate on the case 
P+ = 0, and will not be concerned by the phenomenon 
of eigenvalue crossing. Dropping the subscript for P, we 
consider the Hamiltonian 



P 



1 - IM2 



(61) 



The corresponding equations of motion are 
X P 



N^J2Jl - \Mo 



sin 2 (6lj + V') (62) 



which have to be compared with Eq. ([19[): this is still a 
pendulum equation, but the potential amplitude, depend- 
ing now self consistently on the particles motion through 
M2, explains the presence of untrapped particles, shown 
in Fig. (^). In addition, we have now the right expression 
for the frequency of the fast oscillations lj = w-, and a 
new conserved quantity has been identified: the adiabatic 
invariant P. Using Eq. (p8|), we obtain the A_(M2) rela- 
tion which, together with lo{M2), is perfectly verified by 
numerical simulation as attested by Fig. 

The full efficiency of this procedure is that it pre- 
serves completely the Hamiltonian structure of the prob- 
lem, making it well suited for a statistical treatment, pre- 
sented below. However, it is important to emphasize be- 



Fig. 10. Comparison between the numerical (triangles) and 
analytical (solid line for lj_ given by Eq. (^s])) frequency of 
oscillations. Comparison between the numerical (circles) and 
analytical (solid line for Eq. (^^) amplitude e of oscillation . 
e = 3.610"^ and t„„, = 1500. 



fore, that all fast variables having disappeared, this pro- 
cedure induces a huge gain (of order 1/e) for numerical 
simulations. This allows to study accurately the statistics 
and the dynamics of this effective Hamiltonian up to ex- 
tremely long times, giving conclusion directly relevant to 
the original model. 



4.3 Study of the effective Hamiltonian 

4.3.1 Statistical mechanics of the effective Hamiltonian 

We would like to describe the bicluster as an equilibrium 
state of the effective Hamiltonian (|6l]). Fortunately, due 
to the special form of the potential, which depends only 
on the global variable IM2I, the statistical mechanics of 
the effective Hamiltonian are exactly tractable in the mi- 
crocanonical ensemble. 

There are two conserved quantities: the total energy E 
and the total angular momentum. As the latter only cre- 
ates a global rotation of the system totally decoupled from 
the rest of the dynamics, we can consider the total energy 
without restriction. The density of states Q{E) at a given 
energy, written as an integral in which the energy is split 
in two parts, a kinetic and a potential one, has the follow- 
ing expression 

Q{E) ^ [ dv Qk^niE - V) n^ot{y) 



CC d\M2\f2k,n {E-V{\M2\)) f2,onf{\M2\) (63) 
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Fig. 11. Spacio-temporal evolution of 50 particles. Results 
given by the original Hamiltonian (0) (resp. effective Hamilto- 
nian (pOl)) with (resp. without) the fast small oscillations su- 
perimposed: they are almost indistinguishable. The initial con- 
dition corresponds to particles evenly distributed on the circle 
and with a sinusoi'dally modulated impulsion corresponding to 
an energy density e = 2.5 10~^ and e = 5 10~^. 



where the sign oc accounts for an unimportant constant. 
^kin (resp. f2pot) is the density of states corresponding to 
the kinetic (resp. potential) part of the Hamiltonian, and 
^^conf is the density of angular configuration giving rise 
to a given \M2\: since the potential only depends on IM2I, 
^potiV) is directly proportional to f^conf{\M2\). Their ex- 
pressions are 



^conf{\M2\) 



f^Wdp.sl^-Y^P^-E^ (64) 



N 

Using the classical result for a perfect gas, we obtain 



r{N/2 + 1) 



(66) 



To compute f2confi\M2\), it is possible to use the inverse 
Laplace transform of the Dirac delta function 



5{x) = I dpeP'' 



(67) 



where F is a path in the complex plane running from 
^(p) = to 5(p) = +00. The expression of f2conf 

becomes 



r Jo 



To compute the integral over the angles, we use now the 
Gaussian (also called Hubbard-Stratanovich) transform 



e^P" oc 



+00 



du e *p 



Nf-+Nmu 



This decouples the integration over the angles 9i, and we 
have 



conf{\M2\) cc I dpe 

2-K 



^fl^^^l' / dudve-^^^"+^"^ 



oc I dp e 



-Np\M2 



rdre~^ Io{rf (70) 



where r is the polar radius associated with u and v, and 
Jo is the modified Bessel function of order 0. Using the 
saddle point method on p and r, we have finally 



^ g-A'(2p'|M2|^-ln/o(r*)) 

where p* and r* are implicitly defined by 



*2 

r 

4p*^ 



d ln/o(p) 
dp 



(71) 
(72) 

(73) 
(74) 



Using this result, Eq. (g^ may be evaluated once again 
by the saddle point method, which leads after some sim- 
ple algebra to the equilibrium value IM2I*, defined by the 
following equation 



r*{\M2n = 



1 



4v/l-|M2| 



V2E 



(75) 



-V^-\M2\ 



As shown by Fig. y_2[ this equation has only one solution 
for any ratio E/P; this indicates an absence of phase tran- 
sition. As this solution is non zero, a biclusterization will 
always take place; nevertheless in the large energy limit 
the particles are almost free and IM2I* goes of course to 
zero. 

Once IM2I* is known, it is easy to complete the de- 
scription of the equilibrium state. The temperature T is 
given by 



2 



E- P 



1 - IM2 



(76) 
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Fig. 12. Statistical prediction of the dynamical indicator \M2\ 
as a function of the ratio between the energy E and the adia- 
batic invariant P. 



the distribution of velocities is a Maxwellian with tem- 
perature T, and the distribution of angles has a gibbsian 
shape p(9) cx e~^'^^^/^ with the potential 



v{e) = 



p 



(l-cos(26l + 2V')) (77) 



iV2V2v/l - \M2\- 
This potential is inferred from the equations of motion (|6^) . 

4.3.2 Relaxation to equilibrium 

We have therefore now a complete description of the statis- 
tical equilibrium states of the effective Hamiltonian, gov- 
erned by long range interactions. It has been noticed by 
various authors that the dynamics of such systems may 
sustain long lived metastable states before relaxing to 
equilibrium |l|,|,|5[|2^. Before comparing the "effective 
equilibrium" with the structure created by the dynamics 
of the real Hamiltonian, it is thus necessary to study the 
relaxation to equilibrium. As the relaxation properties of 
long range interacting systems is in itself an important 
problem, we will consider them now. We will in addition 
illustrate which statistical properties of the equilibrium 
distribution is expected to be observed in the real dynam- 
ics, for large N and on time scale reasonable for numerical 
computations. 

Fig. |l^ illustrates the approach to equilibrium of the 
effective Hamiltonian, for three different particle numbers, 
with initially immobile and homogeneously distributed par- 
ticles (this corresponds to the typical initial condition used 



Fig. 13. Relaxation of the first three even moments < M2n > 
(time averaged) toward equilibrium according to numerical 
simulations of the effective Hamiltonian (|6l|. The solid (resp. 
dotted and dashed) lines corresponds to results for A'' = 200 
(resp. iV = 800 and N = 3200) particles. < M2 >, < M4 > and 
< A/e > are represented from top to bottom (< M2n > is a de- 
creasing function of n). They converge toward the equilibrium 
values M| = 0.510, Ml = 0.144 and M^* = 0.028. 



in the real Hamiltonian in |lj^ for instance). We have 
represented < M2„ > (r), the time averages of M2„ from 
initial time to time t, for n equal to 2, 4 and 6. Tem- 
poral fluctuations are thus not visible. We first observe 
finite size effects concerning the equilibrium values : for 
instance Mq converges for the small system N = 200 to a 
value larger than Mg , the equilibrium value. 

More interestingly, we observe also that whereas IM2I 
quickly reaches its equilibrium value, the relaxation time 
of the successive moments IM4I, IMg] strongly depends on 
the system size, and presumably diverges with N . The ac- 
tual dependence of the relaxation time of each moment, 
on the particle number, will not be studied in this paper. 
It would for instance address the issue of whether non 
equilibrium distributions may be observed in the thermo- 
dynamic limit. Such phenomena have indeed already been 
observed in other long range interacting systems p5|,p7t. 

We conclude that for a large, but finite, particle num- 
ber, the moments of the distribution converge toward their 
equilibrium values. When TV increases, the relaxation time 
for the high order moments may however be larger than 
computationally achievable times and, we thus expect the 
simulation to exhibit very slowly evolving non equilibrium 
structures of the effective Hamiltonian. Authors of |4| have 
indeed found a distribution of particles different from the 
Gibbsian shaped distribution ([77|) predicted by the equi- 
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Fig. 14. Comparison of the angular density distribution p{9) 
of particles obtained with the original Hamiltonian (|^) (rep- 
resented with circles) and with the effective one (^l|) (repre- 
sented with the dashed line). Both results have been obtained 
for A'^ = 10'^ particles and are averaged on intermediate times, 



corresponding to r = 10 
Hamiltonian is e = 2.5 10" 



10 . The energy of the original 



librium thermodynamics of the effective Hamiltonian (see 
Fig. I). 

On the contrary our analysis shows that IM2I quickly 
reaches its equilibrium value IM2I*. For this reason, given 
the computational time achievable, in the next section, we 
only use this statistical equilibrium indicator for the study 
of the structure. This will lead to the correct prediction of 
the equilibrium repartition of potential and kinetic energy, 
and of its dependance from initial condition. 



4.4 Comparison with the full Hamiltonian 

If both time scales are clearly separated. Fig. ^ empha- 
sizes the striking agreement between the effective and the 
real dynamics, for short times. We will show in this sec- 
tion that the effective Hamiltonian also provides a good 
description for the long time dynamics. 

For this purpose, we will first compare the evolution of 
the non-equilibrium angular density distribution p(0) for 
both dynamics. The results are reported on Fig. n4 were 
circles show the repartition given by the original Hamil- 
tonian (|l|) whereas the dashed line shows the results of 
the effective dynamics (|6l|). This picture clearly empha- 
sizes that the non equilibrium character of the dynamics 
is fully described by the effective dynamics, and we can 
easily conclude that the effective Hamiltonian reflects very 
well the dynamics of the real Hamiltonian. 



As discussed in the previous section, according to the 
computational times used, only the second momentum of 
the distribution M2 should correspond to the equilibrium 
one. Let us therefore analyze the dependence of the equi- 
librium value of M2 on the initial conditions. The typically 
used initial condition for the numerical simulations in the 
original system (|l|) is an initially quasi homogeneous dis- 
tribution of particles with zero velocity Q . It corresponds 
to the ratio e/^/2P = 1 for the effective Hamiltonian. 
According to Fig. ^ , this implies an equilibrium value 
|Af2|* — 0.51 in perfect agreement with the numerical re- 
sult reported earlier . 

From these results, it is also possible to explain the 
caloric curve T ~ 1.3 e, reported earlier [Q. The total 
energy of the full system is divided in three parts: the 
potential energy of the small oscillations, the kinetic en- 
ergy of the small oscillations and the kinetic energy of the 
slow movements. The two first parts are equal in average 
and form the potential part of the effective Hamiltonian, 
whereas the last one is the kinetic energy of the effective 
Hamiltonian. Using this remark together with the values 
of e/P and |il/2|* at equilibrium, it is easy to derive the 
temperature/energy relationship, using (Uw. In the case 
investigated in Q], we find precisely T = 1.3 e. 

We thus conclude that the dynamics of the effective 
Hamiltonian parametrizes very well the dynamics of the 
real Hamiltonian, for short as well as for long time. This 
allows us to predict statistical properties of the initial sys- 
tem, as for instance the asymptotic value of M2 or the 
partitioning between kinetic and potential energy. More- 
over, let us recall that the effective Hamiltonian gives the 
opportunity to study numerically the relaxation towards 
equilibrium of the bicluster, whereas it was not possible 
in the real dynamics, because of computational limitations 
(let us note that the ratio of the typical time scale of the 
two dynamics is of order 100 or larger). 

All these findings are great successes of this approach. 
Let us nevertheless comment some points that we have 
not yet addressed. 

The first one concerns the limit of validity of our mul- 
tiple time scale analysis. It should be noticed that not 
all initial conditions with small energy would lead to the 
formation of the bicluster. From our analysis, we can con- 
clude that the class of initial conditions leading to this 
formation is the one compatible with the ansatz used: 
however a precise description of this class is not known. 
In particular, we have not study the threshold value of e 
beyond which such a description breaks down. For s going 
toward one, the ansatz ml) will first lead to a nonlinear 
evolution of the slow variable fi in place of the linear one 
given by Eq. (^6|). We think that the critical value Ec above 
which this nonlinear dynamics do not have any more peri- 
odic solution should correspond to the critical value above 
which the bicluster can not form, explaining the transition 
observed in Q. Such an hypothesis should be tested. 

A second point is t he phenomenon of level crossing 
discussed in section 4.2, leading to an exchange of excita- 
tion of the two modes of the system. This phenomenon is 
very similar to level crossing in quantum mechanics adi- 
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abatic problems. It represents a resonance, local in time, 
in which an interaction between the various modes may 
occur, leading to a modification of the adiabatic invariant. 
For sake of simplicity, we have concentrated our work on 
the case of excitation of only the smallest frequency, so 
that such level crossing could not occur. A more general 
study would however be of interest. 

The last point concerns the validity of our approach 
for infinite times. As noted in the introduction, a recent 
paper has shown a long-time degradation of the biclus- 
ter, for a very small number of particles, suggesting its 
transient non-equilibrium nature in that case. For a rea- 
sonable number of particles, our numerical results show 
that such a degradation does not exist, for computation- 
ally achievable times. This point is thus not addressable 
numerically. From a theoretical point of view, results on 
much more simple systems, show that adiabatic invariants 
as described here, should be conserved for a very long time 
(exponential when e goes to 0). Moreover, such very long 
time stability results, would not give any hint on what 
happens for larger times (stability or instability). 



5 Conclusion 

The surprising formation and stabilization of the bi- 
cluster in the HMF dynamics, in contradiction with sta- 
tistical mechanics prediction, is now understood: the small 
collective oscillations of the bulk of particles create an ef- 
fective double- well potential, in which the particles evolve. 

On a first stage the dynamics can be described in 
the context of a forced Burgers' equation. The particles 
trajectories caustics form infinite density regions explain- 
ing an initially diverging density. Early times dynamics is 
therefore very similar to the structure formation in Eu- 
lerian coordinates for a one dimensional self-gravitating 
system ||2^. As in this case, because of the confining po- 
tential, the particles can not move apart as easily after 
the first caustic has been formed, as in the case of the free 
motion on the plane ||2^ . 

On much larger timescales, in order to parametrize the 
fast oscillations of the bulk particles, we have performed a 
variational multiscale analysis. We then have obtained an 
Hamiltonian effective model describing the particle slow 
motion. This description is in very good agreement with 
the initial dynamics. The effective dynamics allows numer- 
ical simulations on time scale much larger than the initial 
dynamics, as the rapid oscillations have been filtered out. 
We have performed the statistical mechanics of this new 
system. The results thus give a statistical explanation of 
the bicluster formation and stabilization. The equilibrium 
properties can then explain the mean value of the second 
moment of the distribution and the repartition between 
potential and kinetic energy. Even if this analysis has only 
been sketched in this work, the effective dynamics is also 
a powerful tool to study the very slow relaxation process 
versus the real equilibrium structure. 

Many physical systems share the main properties of 
this HMF dynamics : very fast oscillations self interacting 
with a slower motion. In addition to the plasma problem 



already cited and always in the context of long range in- 
teracting systems, we may for instance cite the problem of 
interaction of fast inertia gravity waves with the vortical 
motion, for the rotating Shallow Water or the primitive 
equation dynamics, in the limit of a small Rossby Num- 
ber |29) . The main interest of our study is to provide a toy 
model in which such a complex dynamics can be treated 
and analyzed extensively using powerful theoretical tools. 
We point out in particular the usefulness of a variational 
approach in the procedure of averaging the rapid oscil- 
lations. This toy model also permits a clear view of the 
usual problems of such dynamics 

(i) Resonances (or level crossing) in averaging procedures 

(ii) Unusual relaxation processes due to the long range 
nature of the interactions. 

The HMF model is to our knowledge the simplest one in 
which such phenomena occur. The study of these phenom- 
ena is thus a natural extension of our work and should be 
of interest in analyzing more complex systems. 
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Appendix A: Variationnal description of a sim- 
ple adiabatic dynamics 

Let us consider a slowly modulated harmonic oscillator with 
Lagrangian 



L{e,9,t,T) = 



(78) 



2 ' ' 2 

where r — et. We consider the following ansatz 

e{t,T)^ A{T)sinip{t) (79) 

with if = 0(1) and A = edA/dr, but, contrary to the usual 
asymptotic expansion on the equation of motion, we will con- 
sider the variational approach proposed by Witham [^. The 
idea is to separate the two different time scales of the motion 
at the level of the action. The average of the fast variable t 
yields the following effective Lagrangian: 



(80) 



It is thus straightforward to obtain the equation of motion of 
the effective Lagrangian £{A,A,ip,ip). They read 

d_ /d£\ _ dC 
dt [dAJ ^ dA 
d_ (dC\ _ dL_ 

dt \dip j dip 



(81) 
(82) 



= A(<^'-a;^) 
d ( A^if 



dt 



= 



ip = U) 



A (fi = cste 



(83) 
(84) 
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This method emphasizes a new constant of the motion A^ip, 
called adiabatic invariant, that is not exhibited by the usual 
asymptotic expansion on the equations of motion correspond- 
ing to the original Lagrangian (po|). This may be of great in- 
terest in complex systems. Moreover, this invariant necessarily 
appears as the angle </> associated to the fast variable t is always 
cyclic after averaging. The same reasoning also apply for larger 
order ansatzs, showing that an adiabatic invariant should exist 
at any order in e. 

In addition, this method preserves the Hamiltonian charac- 
ter of the problem allowing for instance a statistical mechanics 
treatment as presented in this paper. 

Appendix B: Eigenvalues and eigenvectors of 
the matrix T 

To solve the linear second order equation (^), we need to 
diagonalize the A'^th-order circulant matrix T defined by Tij = 
cos(0i — 9j)/N . We introduce the two vectors X(i/)) and Y(7/j), 
with coordinates Xi{')p) = cos{9i + tp) and Yi{tp) = sm{di + 
tp), and tp an arbitrary phase. T is therefore the sum of the 
projections along these vectors 

T = 1 (X(V)*XW + Y(^)*Y(^)) . (85) 

This proves that the image of T is of dimension two, and that 
only two eigenvalues X± are nonzero. Let us restrict now to the 
two dimensional problem in the Vect(X.,Y) plane (this plane 
does not depend on of course) , and let us choose tp such that 
X and Y are orthogonal, i.e. 

(X.Y) = ^ cos(6li + ip) sm{e, +'tp)=0 . (86) 

i 

This definition of ip leads to 

I A/2 1 = ^ ^cos(26», + 2i>) . (87) 

i 

Once we have the two eigenvectors X+ = X('(/)) and X_ = 
Y (■(/)), their associated eigenvalues are defined by the following 
relationship, valid Vi: 

J2^co8{e^^9,)X±, = X±X±^. (88) 
j 

leading directly to 

= ^ (891 
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